	subroutine viabilidade(PARG,Xcont,Ycont,numpar,numpc,
     &				nfix,viavel,xl_el)

	include 'param_dim.inc'
	logical viavel
	logical intersect
	double precision Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,area
	double precision Xcont(NX),Ycont(NX)
	double precision xl_el,u,v

	double precision PARG(numpc*3),CINC(ninc_max),C0
        double precision x(nparg_max),y(nparg_max),xc,yc,areaT
	double precision Ju,Jv,Juv,compr
        double precision areaD,areaD2 
	integer NNINC(ninc_max)
	integer NUMPI(ninc_max), NUMPCI(ninc_max)
	integer VARIA(ninc_max)


	COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA
	save iteste,areaD
	data iteste/0/


	iteste = iteste +1
	viavel = .true.


! Testa se a orientacao dos pontos de controle esta no sentido horario:
	if(iteste.eq. 1)then
		ipcacum = 0
		do inc=1,NINCL
			do i=1,NUMPCI(inc)
				x(i) = PARG((i+ipcacum)*3-2)
				y(i) = PARG((i+ipcacum)*3-1)
			enddo
		
! 			call centroide(x,y,NUMPCI(inc),0,xc,yc,areaT)

      			call propgeom(x,y,NUMPCI(inc),0,compr,areaT,
     &					xc,yc,Ju,Jv,Juv)
			areaT = - areaT


			ipcacum = ipcacum + NUMPCI(inc)
		
			if(areaT .gt. 0.0)then
				write(*,*)'Erro na orientacao dos pontos de controle'
				viavel = .false.
				return
			endif
		enddo
	endif





! Verifica se o contorno externo estah orientado corretamente
	if(iteste .eq. 1)then
! 		call centroide(Xcont,Ycont,nfix,0,xc,yc,areaD)

      		call propgeom(Xcont,Ycont,nfix,0,compr,areaD,
     &					xc,yc,Ju,Jv,Juv)
		areaD = - areaD

		if(areaD .lt. 0.)then !
			write(*,*)'Malha do contorno orientada incorretamente'
			stop
		endif
		areaD = abs(areaD)
	endif


! Verifica se todos os ptos de controle pertencem ao domínio
	do j=1,numpc
		areaD2 = 0.

		do i=1,nfix
			Ax = Xcont(i)
			Ay = Ycont(i)
			Bx = Xcont(i+1)
			By = Ycont(i+1)
			Cx = PARG(j*3-2)
			Cy = PARG(j*3-1)
			if(i.eq.nfix)then
				Bx = Xcont(1)
				By = Ycont(1)
			endif

			if(area(Ax,Ay,Bx,By,Cx,Cy).lt.0.)then
				write(*,*)'Ponto de controle fora do domínio!'
				viavel = .false.
				return
			endif
		enddo

	enddo




!! Testa se os segmentos de reta formados pelos pontos de controle nao se interceptam 

	ilatacum = 0
	do inc=1,NINCL
		
		ilat1 = ilatacum + 1
		ilat2 = ilatacum + NUMPCI(inc) 
		
		do ilat = ilat1 , ilat2
		
			ia = ilat
			if(ilat.eq.ilat2)then
				ib=ilat1
			else
				ib=ia+1
			endif
	
			Ax = PARG(ia*3-2)
			Ay = PARG(ia*3-1)
			Bx = PARG(ib*3-2)
			By = PARG(ib*3-1)
	! 	    

			jlatacum = 0
			do jinc = 1,NINCL 
				jlat1 = jlatacum + 1
				jlat2 = jlatacum + NUMPCI(jinc)
		    
				do jlat = jlat1,jlat2
					ic = jlat
					if(jlat.eq.jlat2)then
						id=jlat1
					else
						id=ic+1
					endif
			
					Cx = PARG(ic*3-2)
					Cy = PARG(ic*3-1)
					Dx = PARG(id*3-2)
					Dy = PARG(id*3-1)

					if(jlat .ge. ilat+2 .and. id .ne. ia)then

						if(intersect(Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,u,v) .eqv. .true.)then
							viavel = .false.
							write(*,*)'1-INTERSECAO DOS SEGMENTOS ',ilat,jlat
							return
						elseif(u .le. xl_el .or. v .le. xl_el)then
							viavel = .false.
							write(*,*)'PC mto proximos',ilat,jlat
							return
						endif
					elseif(ib .eq. ic .and. jlat .gt. ilat)then

						if(theta(Ax,Ay,Bx,By,Dx,Dy,u,v) .lt. 3.14/18.)then ! 10graus entre as laterais
! 							viavel = .false.
! 							write(*,*)'Angulo theta<10graus ',ilat,jlat
! 							return
						elseif(u .le. xl_el .or. v .le. xl_el)then
							viavel = .false.
							write(*,*)'PC mto proximos',ilat,jlat
							return
						endif
					elseif(ia .eq. id .and. jlat .gt. ilat)then

						if(theta(Cx,Cy,Dx,Dy,Bx,By,u,v) .lt. 3.14/18.)then ! 10graus entre as
! 							viavel = .false.
! 							write(*,*)'Angulo theta<10graus ',ilat,jlat
! 							return
						elseif(u .le. xl_el .or. v .le. xl_el)then
							viavel = .false.
							write(*,*)'PC mto proximos',ilat,jlat
							return
						endif
					endif

				enddo
		    
				jlatacum=jlatacum + NUMPCI(jinc)     
			enddo

		enddo

		ilatacum = ilatacum + NUMPCI(inc)
	enddo
	    
	return 
	end
	

	function theta(Ax,Ay,Bx,By,Cx,Cy,u,v)

	double precision Ax,Ay,Bx,By,Cx,Cy
	double precision theta,u1,u2,v1,v2,u,v

	u1 = Bx-Ax
	u2 = By-Ay
	v1 = Cx-Bx
	v2 = Cy-By
	u = sqrt(u1**2 + u2**2)
	v = sqrt(v1**2 + v2**2)


	theta = acos(abs(u1*v1+u2*v2)/(u*v))
	
! 	write(*,*)'u1,u2,u',u1,u2,u,'v1,v2,v',v1,v2,v,'theta ',
!      &			acos((u1*v1+u2*v2)/(u*v))


	return
	end

	
	function intersect(Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,u,v)
	
	logical intersect,collinear,Xor,left
	double precision Ax,Ay,Bx,By,Cx,Cy,Dx,Dy
	double precision u1,u2,v1,v2,u,v
	
	if (collinear(Ax,Ay,Bx,By,Cx,Cy) .or. 
     &	    collinear(Ax,Ay,Bx,By,Dx,Dy) .or. 
     &      collinear(Cx,Cy,Dx,Dy,Ax,Ay) .or. 
     &	    collinear(Cx,Cy,Dx,Dy,Bx,By))then
	    
	    intersect = .false.
	else
	    intersect= (Xor( left(Ax,Ay,Bx,By,Cx,Cy), 
     &                       left (Ax,Ay,Bx,By,Dx,Dy) ) .and.  
     &			Xor( left(Cx,Cy,Dx,Dy,Ax,Ay), 
     &			     left(Cx,Cy,Dx,Dy,Bx,By) ) )
	endif

	u1 = Bx-Ax
	u2 = By-Ay
	v1 = Dx-Cx
	v2 = Dy-Cy
	u = sqrt(u1**2 + u2**2)
	v = sqrt(v1**2 + v2**2)

	
	return
	end
	
	
	function collinear(Ax,Ay,Bx,By,Cx,Cy)

	logical collinear
	double precision Ax,Ay,Bx,By,Cx,Cy,area
		
	if( area(Ax,Ay,Bx,By,Cx,Cy) .eq. 0.)then
	    collinear = .true.
	else
	    collinear = .false.
	endif

	return
	end

	function area(Ax,Ay,Bx,By,Cx,Cy)

	
	double precision Ax,Ay,Bx,By,Cx,Cy
	double precision area
	double precision area_parc

	area =  area_parc(Ax,Ay,Bx,By) +
     &		area_parc(Bx,By,Cx,Cy) +
     &		area_parc(Cx,Cy,Ax,Ay)



! 	area = ( ax * by - ay * bx +
!      &           ay * cx - ax * cy +
!      &           bx * cy - cx * by ) !retorna area*2
     
!      	write(*,*)'areas:', area1,area2
	return
	end
	
	function left(Ax,Ay,Bx,By,Cx,Cy)
	
	double precision Ax,Ay,Bx,By,Cx,Cy,area
	logical left
	
	if(area(Ax,Ay,Bx,By,Cx,Cy) .gt. 0.)then
	    left = .true.
	else
	    left = .false.
	endif
     
	return
	end
	
! 	logical function Xor(x,y)
! 	
! 	logical x,y
! 	
! 	if((x .eqv. .false. .and. y .eqv. .true.) .or. 
!      &	   (x .eqv. .true. .and. y .eqv. .false.))then
! 	    Xor = .true.
! 	else
! 	    Xor = .false.
! 	endif
! 	
! 	return
! 	end
	    
